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Abstract -Bilayer lipid membranes (BLMs) are an essential component of many biological sys¬ 
tems, forming a functional barrier between the cell and the surrounding environment. When 
the membrane relaxes from a structural perturbation, the dynamics of the relaxation depends on 
the bilayer structure. We present a model of a BLM in a viscous solvent, including an explicit 
description of a Thick’ membrane, where the fluctuations in the thickness of a monolayer leaflet 
are coupled to changes in the lipid density within that monolayer. We find dispersion relations 
describing three intuitive forms of bilayer motion, including a mode describing motion of the in¬ 
termonolayer surface not noted previously in the literature. Two intrinsic length scales emerge 
that help characterise the dynamics; the well known Saffman-Delbriick length and another, C, 
resulting from the intermonolayer friction. The framework also allows for asymmetry in the BLM 
parameters between the monolayer leaflets, which is found to couple dynamic modes of bilayer 
motion. 


Introduction. — Bilayer lipid membranes (BLMs) are 
ubiquitous in nature and common in practical applica¬ 
tions. The complexity of biological membranes is typified 
by the cell membrane, a complex dynamic mix of pro¬ 
teins and lipids that act as the gatekeeper to the cell, the 
lipids additionally providing an active host for membrane 
proteins [T]. The double layer structure of a BLM allows 
for complex equilibrium conformations and rich dynamical 
behaviour [2]. The BLM’s role as a boundary between two 
diverse environments means that the membrane composi¬ 
tion will rarely be symmetric across the monolayer leaflets 
that comprise the bilayer; the two membrane surfaces may 
have different functional requirements. An asymmetry 
in composition will lead to an asymmetry in behaviour, 
which can be used as a control mechanism in biological 
processes Hi]. The symmetry of the membrane can also 
be broken by the imposition of an external field, such as 
the electric fields used in electroporation [J. The ways in 
which asymmetry may effect and even assist the biological 
functions of a BLM is as yet unclear. Equally as unclear 
is the complex dynamic behaviour that can emerge in an 


asymmetric BLM. 

At the simplest level, a BLM can be modelled math¬ 
ematically as a thin sheet embedded in fluid [a El , as 
employed by Brochard & Lennon to model the fluctua¬ 
tions of red blood cell membranes HO]. A thin sheet can 
only generate a single hydrodynamic mode, where the dis¬ 
placement of the membrane is resisted by the curvature 
and tension of the membrane and energy is dissipated in 
the surrounding fluid. A similar approach of coupled thin 
sheets has been used by Lu and Cates to model the hy¬ 
drodynamics of surfactant films m- 

Seifert & Langer [Ellis] expanded the description to 
a BLM comprising two compressible monolayer leaflets 
coupled to the bending of the bilayer. The movement of 
lipids within a monolayer is restricted by intermonolayer 
friction, the viscosity of the surrounding solvent and the 
viscosity within a monolayer. Seifert & Langer’s model 
considers three modes of motion; a height mode (as pre¬ 
viously modelled by Brochard & Lennon uni) , an aver¬ 
age density mode, where the monolayer densities move in 
phase, and a density difference mode, where the mono¬ 
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layer densities move out of phase. The average density 
mode is not considered within the dissipative dynamic 
framework and is considered a purely propagating mode, 
which leaves combined dispersion relations for the re¬ 
maining dynamic modes, where the dominant relaxation 
mode depends on the length scale of the bilayer move¬ 
ment. This model shows quantitative agreement with ex¬ 
perimental work UMI] and has subsequently become the 
most widely used framework for studies of BLM dynamics. 
Hdmberg & Muller [18] generalised the Seifert & Langer 
theory to include the bilayer’s inertia and surface ten¬ 
sion although the inclusion of inertial effects does not im¬ 
prove the agreement with experimental results since most 
nanoscale BLM undulations occur at vanishing Reynolds 
number. Coarse-grained simulation models of membranes 
allows for the direct calculation of fluctuation spectra, 
where the Seifert & Langer framework has again been used 
to effectively describe the fluctuations observed 1121120]. 

Recent computational EDna and experimental [231125] 
work has highlighted the role fluctuations in the membrane 
thickness can play in the dynamics. Thickness fluctua¬ 
tions are not explicitly considered in the Seifert & Langer 
model, as the BLM is physically represented as a thin 
sheet. Fluctuations in the lipid density in each leaflet are 
analogous to thickness fluctuations, but neglect the move¬ 
ment of fluid due to thickness changes, which will affect 
the dynamics. Advances in computer simulation have also 
led to debate around the Seifert & Langer model of BLM 
tension, suggesting that BLM tension can arise from sev¬ 
eral sources and can even vary between monolayer leaflets 
in apparently Tensionless’ BLMs [26]. While this has led 
to new BLM models based on renormalized elastic theory 
m or a viscoelastic tensor-based approach [28] , these are 
not conclusive and the issue is still unresolved. 

Here we present a new model of membrane dynamics 
that includes a physical description of a ‘thick’ membrane 
comprising two monolayers, with the freedom to vary pa¬ 
rameters between the monolayers. The thickness is ex¬ 
plicitly included by evaluating the boundary conditions 
at the membrane surface, rather than on the membrane 
midplane, as previously considered. We describe the dy¬ 
namics in terms of intuitive modes of motion of the bilayer, 
including a form of bilayer relaxation not noted in previ¬ 
ous models. For a symmetric bilayer, the model generates 
dispersion relations that can be studied analytically and 
are associated with the modes of bilayer motion. 

Method. — 

Geometry. We model an asymmetric BLM by vary¬ 
ing the monolayer thickness, membrane viscosity and area 
compressibility between the monolayer leaflets, which cou¬ 
ples the modes of bilayer motion, complicating the dissi¬ 
pative relaxation of an asymmetric BLM. We consider a 
planar BLM suspended in viscous fluid, with the unper¬ 
turbed bilayer normal parallel to the 2 ;-axis and intermono¬ 
layer surface spanning the x-y plane at z = 0 , as shown in 
Fig. m The upper (-f) and lower (—) bilayer surfaces are 


described by the height functions (r) and h- (r): 

h+ (r) = d+ (r) + s (r) (la) 

/I-(r) =-d_ (r) + s (r), (lb) 

where d± (r) is the thickness (strictly positive) of a mono- 
layer and s (r) is the height of the intermonolayer surface. 
We define r = (x^y). This explicit description of a ‘thick’ 
bilayer is a distinctive feature of our model and better 
reflects the physical nature of a BLM. 



Fig. 1: A 2D projection of the schematic geometry of our bi¬ 
layer model. 

Free Energy. The BLM description contains three de¬ 
grees of freedom; the two monolayer thicknesses and the 
height of the intermonolayer surface. By comparison, the 
Seifert & Langer model considers two degrees of freedom, 
the bilayer height and out of phase variation of the mono- 
layer density, where the average density mode is neglected 
within the non-inertial framework [12] . The free energy of 
each monolayer F± is modelled after the Helfrich-Canham 
Hamiltonian in the Monge representation [29] with, in 
principle, different elastic constants for each monolayer: 


F — Fj^ -|- F- + F frame] 


where. 


F± = 


= \ j [Kb± ( v 2 / i ±)"+ 7 ,± {Vh±f 


-\-nA± 


d± - do± 
dozL 


do± vs 


dA, 


Fframe = 7/. [V(/i+ + /i_)]V2, 


(2a) 


(2b) 

(2c) 


where is the membrane bending rigidity, do± is the 
unperturbed monolayer thickness and ka is the area com¬ 
pressibility. The second term in the area compressibility 
couples the curvature of the membrane and the thickness 
of the monolayers m- The surface tension 7 ^ restricts 
variations in the monolayer/water interfacial area while 
we introduce an additional frame tension term to restrict 
changes in the total membrane area. The contribution 
to the free energy from the frame tension jfr cannot be 
split between leaflets and is a property of the whole bi¬ 
layer. This description of the BLM covers every form of 
deformation of the bilayer. 
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Dynamics. The surrounding fluid is assumed to be 
incompressible (V • v = 0) and in the non-inertial Stokes 
limit (as used by Seifert & Langer m), where v is the fluid 
velocity [18]. This simplifies the corresponding Navier- 
Stokes equation to 


The final boundary condition ensures continuity of the 
normal component of the velocity across the monolayer- 
monolayer interface 

dth+- dtd+= dth-- dtd-. (9) 


Vp = r]V‘^\. (3) 

where p is the fluid pressure and rj is the viscosity. Each 
monolayer obeys the continuity equation, 

^ + V_L ■ (nvi) = 0 (4) 

where is the (two-dimensional) velocity of the lipids and 
n is the two dimensional number density of lipids within 
the monolayer. The corresponding area per lipid is a = 
1/n. We can assume that the volume vq = ad per lipid 
is fixed as this will relax much faster than perturbations 
to the area or thickness [27] • For small perturbations in 
a,d, and v/, the continuity equation, linking monolayer 
thickness changes with the flow of lipids within the leaflet 
becomes; 

dtSd± = -do±V±_ ■ vf (5) 

to lowest order in all perturbations, where do is the un¬ 
perturbed monolayer thickness and Sd is the perturbation 
to that thickness. 

Boundary Conditions. The surrounding fluid exerts a 
shear stress on the monolayer-fluid interface, denoted 
where i = x^y. Friction also arises at the inter monolayer 
surface due to the relative velocities of the two monolay¬ 
ers. We treat the BLM as coupled two-dimensional fluid 
monolayers, each of which obeys a Navier-Stokes equa¬ 
tion. Including body forces and assuming Stokes flow, the 
Navier-Stokes equation for each monolayer is 

- do±Vi^ + ± T± T b (v+ - v",) = 0, (6) 

where b is the coefficient of intermonolayer friction and /i± 
is the viscosity of the monolayer. The first term represents 
the pressure within each monolayer. In the Stokes approx¬ 
imation the forces across each interface (monolayer/fluid 
and monolayer/monolayer) balance, which leads to the 
normal force balance; 


The dynamics are fully specified by Navier-Stokes equa¬ 
tions and the boundary conditions and we require that the 
velocity and any deviation in the pressure in the surround¬ 
ing fluid vanish far from the bilayer (v^ ^ 0 as z ^ ± 00 ). 

Calculation of Dynamic Modes. We next recast the 
bilayer description into collective degrees of freedom that, 
together with the intermonolayer surface height 5 , reflect 
all possible modes of bilayer movement: 

h = -t- (10) 

u = — h-. 


Here u describes peristaltic undulations in the bilayer and 
h describes whole bilayer undulations. If h and s undulate 
with the same wavelength and amplitude, then the bilayer 
undergoes a thickness-preserving undulation, comparable 
to the height undulations of previous dynamic models m 
[ 12 ]. This recasts the free energy framework in terms of h, 
u and 5 , rather than d±/h± and s. 

We expand in a Fourier series, 

= ( 11 ) 


and similarly for u (r) and s (r). Here q = (qx^Qy) is the 
Fourier wavelength of the undulation. The dynamical vari¬ 
ables then obey 



where 



Mil 0 0 

M21 M22 M23 

M31 M32 M33 


( 12 ) 


(13) 


The matrix Mq contains all the fundamental information 
of the system (to linear order). The solutions are given by 


Ai{t)=Aoeie where i = 1, 2, 3. (14) 


6F SF , 

-^-= T+ 

Sh^ Sh- 


(7) 


We prescribe non-slip and impermeable boundary condi¬ 
tions at each monolayer-fluid interface at 2 : = Fho. The 
boundary conditions are evaluated at the membrane sur¬ 
face, in contrast to previous studies, where they are eval¬ 
uated at membrane midplane [l2l[T8] . This more physical 
description will capture the membrane thickness: 


vf (r, 2 ; = ±ho) = ±dth± (r) (8a) 

v± (r, z = ±ho) = v± (r). (8b) 


where Aq is the initial amplitude of the bilayer deforma¬ 
tion represented by the ith eigenmode, which has eigen¬ 
value (decay rate or dispersion relation) uji (g), and eigen¬ 
vector e^. 

Results. — Since the matrix Mq is symmetric in qy 
we assume g^ = 0 without loss of generality. We first 
consider a symmetric bilayer, where the BLM parameters 
are identical in each monolayer leaflet. This simplifies Mq 
into diagonal components and one off-diagonal component 
(M 31 ) which couples the whole bilayer mode, h and in¬ 
termonolayer surface height s. The diagonal components 
provide analytic expressions for the eigenvalues uji (g). 
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Table 1: Table of parameters. 


Symbol 

Value 

Reference 

Kb 

40 ksT 

m 

Is 

0.1 mNm“^ 

m 

Ifr 

0.3 mNm“^ 

m 

ka 

0.14 Nm-i 

m 

do 

2 nm 

m 

R 

10 nPa ms 

m 

77 

1 mPas 

m 

b 

10 MPasm“^ 

m 


Peristaltic mode. The first mode has eigenvalue e = 
Uq and is a pure peristaltic dynamic mode describing un¬ 
dulations in the bilayer where the thicknesses of the mono- 
layers undulate in phase. The decay rate is 


{qx) = M 22 = 


Qx I^A 

f^Qx T 2 77 


Qx i^A 

V {^SdQx + 2) ’ 


(15) 


in which the area compressibility is damped by the vis¬ 
cosities of the solvent and lipids. The effect of the mem¬ 
brane viscosity can be understood when the eigenvalue 
is expressed in terms of the Saffman-Delbriick length 
^SD = t^/'Q [S]. The Saffman-Delbriick length repre¬ 
sents the length scale below which the membrane viscos¬ 
ity dominates the surrounding fluid viscosity. For wave¬ 
lengths larger than IsD {q^sd ^ 1) the fluid drag domi¬ 
nates while for shorter wavelengths the in-plane monolayer 
viscosity dominates and the damping is independent of 
wavelength. For the membrane viscosity (/i ~ lOnPams) 
and the fluid viscosity (77 ~ 1 mPa s) that are used here 
isD ~ 0.1/im therefore uJu{q) — ^a !for smaller scale 
undulations. The Saffman-Delbriick length emerges from 
considering the surrounding fluid moved by thickness fluc¬ 
tuations, an effect neglected previously nans]. The limit¬ 
ing value of this dispersion relation (at high wavenumber 
{q) ^ 10^s“^) is comparable to the experimental relax¬ 
ation rate (100ns) for thickness fluctuations measured by 
Woodka et al [23] . 

This dispersion relation is comparable to the density 
difference region of an eigenvalue found by Seifert and 
Langer [T2j[T3]. Seifert and Langer found a mixed eigen¬ 
value where terms representing the peristaltic bilayer dy¬ 
namic mode dominate for a particular region of q values 
given in their paper as 72 ioi q^ 2 tika! {b {i^b + 2 ^q/^a)), 
as the combined density and height variables used are cou¬ 
pled when describing this form of bilayer motion. We find 
a single eigenvalue due to the decomposition of membrane 
motion into the collective degrees of freedom (Eg. [TO]) . 

Ripple mode. We next consider a mixed dynamic 
mode that couples the bilayer mode h and the internal 
surface mode s, although the contribution from the lat¬ 
ter has a very weak wavelength dependence and is slightly 


suppressed at low q. The eigenvector is given by 


^hs — 


1 

0 

Mil—M33 - 


where 


M31 


Mss 


Mn 

2 


(16) 


In this dynamic mode, the internal surface s undulates 
in phase with the bilayer surface in order to preserve the 
thickness of each monolayer, and the decay rate is 


^hs {<lx) = Mil 


Cis + 27/r + Kbql)qxe 
271 


(17) 


Relaxation is driven by the membrane bending rigidity 
(Kb) and the surface and frame tensions, and resisted by 
the fluid viscosity. This relation is equivalent to the ripple- 
type dynamic mode of Brochard & Lennon m- Here the 
dynamic mode is embodied in a single eigenvalue, whereas 
in previous models the dynamic mode dominates within a 
larger, more complex eigenvalue (the ’fast’ mode in Seifert 
& Langer) in the shorter wavelength range q ~ 10 “^ — 
lnm“^. 

Internal mode. Finally, we find a dispersion relation 
associated with movement of the internal surface Sq] 


ojs (qx) = M33 


ka ql {ql (Tq - 1) 
ri{2qx + qx^sD + 2 £sd/£1) ’ 


(18) 


where the length scale £r = ijnfb, compares the inter- 
monolayer friction and the membrane viscosity. This dis¬ 
persion relation is a balance between the area compress¬ 
ibility (i^a) and the viscosities 77 and /i and the inter¬ 
monolayer friction, b. A similar dependence is seen in the 
‘slipping’ region (g ^ 10“^ — lnm“^) of the slower mode 
of Seifert & Langer (denoted by 72 in their work). The 
length scale ir represents the size at which the forces from 
the membrane viscosity and the intermonolayer friction 
on a monolayer inclusion balance. For a typical bilayer 
(/i ^ lOnPams, b ^ lOMPasm”^ ) ir ~ lOnm isD as 
the intermonolayer friction is much larger than the mem¬ 
brane viscosity, although current estimates of b do show 
large variation m- The movement of lipids within the 
monolayer also drags the surrounding solvent. The disper¬ 
sion relation tends to zero as gdo ^ 1, where the Fourier 
wavelength approaches the membrane thickness. This un¬ 
physical behaviour is expected in this limit, as the Stokes 
approximation begins to break down. At short wave¬ 
lengths, the inclusion of inertial effects have been found to 
lead to a dependence [32] . The Stokes approximation 
breaks down for only the internal mode uUs as this contains 
the least dissipation of motion in the surrounding fluid of 
all modes and hence the inertia/dissipation balance that 
leads to the Stokes approximation is disrupted before it 
will be in the remaining modes. This relevant portion of 
the figure is marked on figures [2] & |3l 
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Fig. 2: Dispersion relations (uj (q)) as a function of the Fourier 
wavenumber (q). The wavenumbers qsD = and qr — 

2'Klir associated with the Saffman-Delbriick length Isd and 
the new length scale Ir are shown. Power laws are shown 
for comparison. The shaded region of ujs indicates where the 
Stokes approximation breaks down, leading to the false impres¬ 
sion of a maximum value. The parameters used to generate the 
figure are shown in Table [T] 

The dispersion relations are shown in Fig. [2j The dis¬ 
persion relation ujg associated with the internal mono¬ 
layer/monolayer surface has the smallest magnitude for 
the entire range of q^ meaning that this type of bilayer 
motion will be the slowest to relax. At long wavelengths 
the fastest mode is the peristaltic mode uju, before cross¬ 
ing over to the thickness preserving ripple mode ujj^g. This 
crossover occurs at q = 10“^ nm“^ equivalent to an undu¬ 
lation wavelength of 0.1 //m, which is comparable with the 
crossover value between the ‘slipping’ and ‘ripple’ dynamic 
modes found in previous models [lain], simulations [33] 
and experiments [20] . 

The eigenvalues and uju exchange their relative mag¬ 
nitudes as the Fourier wavelength q is increased, but the 
eigenvectors show an extremely weak dependence on q and 
so that the nature of each mode does not change as q in¬ 
creases. In the model first presented by Seifert & Langer 
[T2j[T3], two dynamic modes interchanged their character 
depending on q range, which makes it difficult to interpret 
the bilayer motion. 

Asymmetric BLM. The results considered so far 
are for a symmetric BLM, which is clearly an idealised 
case, since biological membranes are often asymmetric [Ij. 
Varying the bending modulus, Hh and the stretching mod¬ 
uli, 76 and is non-trivial. Here we vary the area 
compressibility, the relative leaflet thickness, do, and 
the membrane viscosity, /i, between leaflets. In this case, 
the matrix Mq does not simplify, and expressions for the 
eigenvalues do not lead themselves to analytical inspec¬ 
tion, with the exception of Mu, which is unchanged. The 


exponential in Mu is unaffected as the do comes from 
a boundary condition which depends on the whole bilayer 
thickness, while we consider variation in the relative leaflet 
thicknesses. The changed terms of the matrix are given 
by 

M2i = qx [ KA - a ‘^ {‘^b + 27]qx +^ i + ql ) ( 19 ) 

-ka+ (2 6 + 27]qx + l^-ql)] /2 Ci 
M 22 = -qx [K,A-a^ (2 6 + 2-qq^ + ji+ql) ( 20 ) 

+ka+ (2 6 + 2r]qx + IJ-ql)] /2 Ci 
M23 = qx [KA-a'^ {dl-ql - l) (2 6 + 2773^, + i^+ql) (21) 
-ka+ {a^ dl_ql - l) (2 6 + 2-qq^ + ^l-ql)] /2 Ci 
M31 = {q^ ri [^A-a^ (l^+qx + 2 v) (22) 

+ka+ {li-qx + 2/?)] 

- a ^ (7s + 2'yfr + Kbql ) C'l} /4 77 Ci 
M32 = ql [a^HA- (2 77 + M+fe) (23) 

-ka+ (2 ?7 + H-qx)] /4 Cl 

^33 = ql [ oP ' KA - { dl_ql - 1 ) (2 77 + ji + q ^) ( 24 ) 

+ KA + { a ^ dl_ql - 1) (2 77 + H - qx )] I^Ci 

where 

a = do +/ do - ( 25 ) 

Cl = [4 677 + ^2, (77-+77+) (2 77^2,+ 6) ( 26 ) 

+4772^^ + 77+77-9^] . 

The numerically calculated eigenvalues are shown in figure 

El 



Wavenumber ^(nm 

Fig. 3: A log-log plot of the dispersion relations (cc (^)), show¬ 
ing the effect of increasing asymmetry. The dispersion relation 
for the ripple mode undulation {q) is unchanged by the 
asymmetry considered in this study. The parameters used to 
generate the figure are shown in Table [U 

The position of the crossover wavelength between 
(q) and uJu {q) depends upon the degree of asymmetry. 
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as do the dispersion relations Uu (q) and ujg (q). Introduc¬ 
ing asymmetry has a seemingly small effect on the mag¬ 
nitude of the dispersion relations, but does inffuence the 
character of the modes, as seen in the eigenvectors shown 
in figure [3l The dispersion relations ujg {q) and uju (q) be¬ 
come mixed modes associated with both u and 5; so that 
a peristaltic excitation of an asymmetric BLM will relax 
both by peristaltic undulations {u) and by movement of 
the internal surface (s). If the dynamic modes are cou¬ 
pled, then an asymmetric membrane would relax at the 
rate of the fastest dynamic mode, as this would provide 
the fastest return to equilibrium. This changes the dissi¬ 
pative dynamics of BLM in experimental membranes. 


Conclusions. — We have constructed a new model 
and framework for the dynamics of BLMs that generalises 
to asymmetric membranes. The framework provides an 
intuitive description of both the BLM and the forms of 
bilayer motion. We explicitly include the effect of thick¬ 
ness fiuctuations on the surrounding fiuid and find a lim¬ 
iting value in the dispersion relation associated with these 
thickness fiuctuations that agrees well with recent experi¬ 
ments [23] . Two length scales emerge from the dispersion 
relations: the previously well studied Saffman-Delbriick 
length isD^ and a new lengthscale, ir a monolayer equiv¬ 
alent to the Saffman-Delbriick length which balances in¬ 
ternal monolayer viscosity with internal friction. The ad¬ 
ditional length scale emerges as the associated form of 
bilayer motion (s) is considered uncoupled from the other 
forms of bilayer motion. Asymmetry does not significantly 
effect the relaxation dynamics of the bilayer but instead 
couples the dynamic modes, meaning an asymmetric mem¬ 
brane (e.g. most biological membranes) will relax by a 
combination of dynamic modes. An asymmetric surface 
tension or bending rigidity may change the dynamics more 
significantly but represents a significant extension to the 
model. Biological membranes contain numerous embed¬ 
ded molecular complexes. The effect of membrane be¬ 
haviour on these complexes has been studied elsewhere 
[3^ . and their inherent asymmetry (either by their con¬ 
formation or by their interactions with the surrounding 
lipids) will affect membrane dynamics. Complex embed¬ 
ded structures such as a transmembrane ion channel could 
be modelled within our existing membrane framework, but 
would be best treated within a fully three dimensional 
description, which would make obtaining analytic solu¬ 
tions difficult. The perfect symmetry required to observe 
pure hydrodynamic modes presents a challenge for exper¬ 
imentalists. However, the change in the wavelength of the 
crossover between the slowest relaxing modes could pro¬ 
vide an opportunity for experimental verification. Recent 
advances in experimental methods allow for the prepara¬ 
tion of intentionally asymmetric membranes Obser¬ 
vations of the fiucuation spectra for these membranes may 
validate the predictions of the model. 
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